Purpose: To measure the growth rate of all strains, starting with cultures grown to mid log phase

Protocol:

library(tidyverse)
library(xlsx)
library(broom)
library(rstatix)
library(ape)
library(formattable)
library(corrplot)
library(ggbeeswarm)
library(cowplot)
library(ggpubr)
library(grid)
library(gridGraphics)

`%!in%` <- negate(`%in%`)

Short cuts and file paths:

strains <- c("14018", "315-A", "C0011E4", "JCP7275", "C0084H9", "JCP8108", "JCP8017B", "C0040C2" , "UM35",  "JCP8066", "C0093B3", "AMD", "C0096A1", "C0102A2",  "UM224", "Gv5-1", "C0179B4", "C0056B5", "Gv101", "C0100B2", "C0101A1", "CMW7778B")

strains0 <- c("14018", "315-A", "C0011E4", "JCP7275", "C0084H9", "JCP8108", "JCP8017B", "C0040C2" , "UM35",  "JCP8066", "C0093B3", "AMD", "C0096A1", "C0102A2",  "UM224", "5-1", "C0179B4", "C0056B5", "101", "C0100B2", "C0101A1", "CMW7778B")

species <- c("Gardnerella vaginalis", "Gardnerella sp. 2", "Gardnerella sp. 3", "Gardnerella piotii", "Gardnerella leopoldii", "Gardnerella swidsinskii", "Gardnerella sp. 7", "Gardnerella sp. 8", "Gardnerella sp. 10", "Gardnerella sp. 11", "Gardnerella sp. 12")

batches <- c("1A1", "1A2", "1B1", "1B2", "2A1", "2A2", "3A1", "3A2")

figureOut <- "../../experiments_figures"
suppTableOut <- "../../experiments_figures/supplementary_posthoc_tables.xlsx"

Import data

midlog_raw_data <- read_csv("./data/gc_midlog_raw_data.csv")
curve_raw_data <- read_csv("./data/gc_curve_raw_data.csv")

Plot settings

theme_set(theme_bw())

speciesColor <- scale_color_manual(values=c("#CC79A7", "#D55E00", "#56B4E9", "#0072B2", "#E69F00", "#009E73", 
                                                 "darkgreen", "gold3", "darkred", "darkblue",
                                                 "mediumpurple4"))

speciesFill <- scale_fill_manual(values=c("#CC79A7", "#D55E00", "#56B4E9", "#0072B2", "#E69F00", "#009E73", 
                                                 "darkgreen", "gold3", "darkred", "darkblue",
                                                 "mediumpurple4"))

Set up growth

Propagate twice and then incubate to mid log. 1:10 Dilution for mid Log: Protocol 1: 4 hour incubation Protocol 2: 7 hour incubation Protocol 3: 15 hour incubation ## Clean data

midlog_df <- midlog_raw_data %>%
  group_by(Number, BioRep, Species, Strain, Batch, Step) %>%
  mutate(OD=OD-mean(.$OD[.$Number=="Blank"]),
                      OD=case_when(OD>0~OD,
                               OD<=0~0)) %>%
  filter(Number!="Blank") %>%
  mutate(Protocol=str_sub(Batch, 1, 1),
         BatchNum=str_sub(Batch, 3,3),
         Step=factor(Step, levels=c("Overnight 1", "Overnight 2", "Dilution Start", "Mid Log")),
         Strain=factor(Strain, levels = strains),
         Species=factor(Species, levels = species)) %>%
  with_groups(c("Number", "BioRep", "Species", "Strain", "Batch", "Step", "Protocol", "BatchNum"), summarize, ODm=mean(OD), ODsd=sd(OD)) %>%
  filter(Step!="Dilution Start") %>%
  mutate(StrainN=as.numeric(Strain),
         rectFill=case_when(StrainN %% 2 == 0 ~ "A",
                            StrainN %% 2 != 0 ~ "B")) %>%
  ungroup %>%
  mutate(Sample=paste(Strain, Batch, BioRep, sep="."))

Assess

Assess mid log growth for removing poor growers with OD600 < 0.05

midlog_df %>%
  filter(Step=="Mid Log") %>%
  ggplot(aes(x=Strain, y=ODm, shape=BioRep, color=Batch)) +
  geom_point() +
  geom_hline(aes(yintercept=0.05), linetype=2, color="gray") +
  theme(axis.text.x = element_text(angle=45, hjust=1)) +
  labs(x=NULL, y=bquote(OD[600]), shape="Biological Replicate", color="Experiment Batch")

Create quality filter for OD600 > 0.05

midLogFail <- midlog_df %>%
  filter(Step=="Mid Log", 
         ODm < 0.05) %>%
  .$Sample
  
midLogFail
##  [1] "C0100B2.1A1.A" "C0100B2.1A2.A" "C0100B2.1A1.B" "C0100B2.1A2.B"
##  [5] "C0096A1.2A2.B" "JCP8066.1A2.A" "JCP8066.1A2.B" "UM224.1A1.B"  
##  [9] "C0100B2.1B1.A" "C0100B2.1B1.B"

Apply quality filter and plot set-up growth endpoints

# main plot with quality filtering
overnightPlot <- midlog_df %>%
  filter(Sample %!in% midLogFail) %>%
  ggplot(aes(x=Strain, y=ODm)) +
  geom_rect(aes(xmin=StrainN-.5, xmax=StrainN+.5, ymin = 0, ymax = 1, fill=rectFill), alpha = 0.2, show.legend = FALSE) +
  geom_pointrange(aes(ymin=(ODm-ODsd), ymax=(ODm+ODsd), color=Species), position = position_quasirandom(), size=0.3, alpha=0.7, show.legend = FALSE) + 
  speciesColor +
  scale_fill_manual(values=c("white", "gray"), labels=c("A","B")) +
  ylim(0,1) +
  facet_grid(.~Step) +
  coord_flip() +
  scale_x_discrete(limits=rev) +
  theme(axis.text.x = element_text(size=9),
        axis.text.y = element_text(size=8.5),
        axis.title = element_text(size=13),
        strip.text = element_text(size=13),
        legend.text = element_text(size=11),
        panel.grid = element_blank(),
        legend.position = "none") +
  labs(x=NULL, y=bquote(OD[600]), shape="Experiment Batch")

#annotation bars
protocolBar <- midlog_df %>%
  select(Strain, Species, Protocol) %>%
  unique %>%
  ggplot() +
    geom_bar(aes(x=Strain, y=2, fill=Protocol), stat="identity", width = 1.1) +
    scale_x_discrete(limits=rev) +
    scale_fill_manual(values=c("cornflowerblue", "burlywood1", "forestgreen")) +
    coord_flip() +
    theme_void() +
    theme(legend.position = "none")

speciesBar <- midlog_df %>%
  select(Strain, Species, Protocol) %>%
  unique %>%
  ggplot() +
    geom_bar(aes(x=Strain, y=1, fill=Species), stat="identity", width=1.1) +
    scale_x_discrete(limits=rev) +
    speciesFill +
    coord_flip() +
    theme_void() +
    theme(legend.position = "none")

protocol_legend <- ggpubr::get_legend(protocolBar +
    guides(color=guide_legend(ncol = 1)) +
    theme(legend.position = "bottom",
          legend.direction = "horizontal"))

species_legend <- ggpubr::get_legend(speciesBar +
    guides(color=guide_legend(ncol = 1)) +
    theme(legend.position = "bottom",
          legend.direction = "horizontal"))

a <- plot_grid(overnightPlot, speciesBar, protocolBar, nrow = 1, align = "h", axis = "bt", rel_widths = c(1, 0.02, 0.02))
b <- ggdraw(plot_grid(a, as_ggplot(species_legend), nrow=2, rel_heights = c(1, 0.3)))
setUpGrowthFigure <- cowplot::plot_grid(b, as_ggplot(protocol_legend), ncol=1, nrow=2, rel_heights = c(1, 0.1))
setUpGrowthFigure

#ggsave(file.path(figureOut, paste(Sys.Date(), "Figure2_curves_setup_growth.png", sep="_")))

Growth Curves

Assess blanks

blanks0 <- curve_raw_data %>%
  filter(Number=="Blank") %>%
  select(well, Number, Time, OD, Batch)
  
blanks0 %>%
  ggplot(aes(x=Time, y=OD, color=well)) +
  geom_point() +
  facet_wrap(~Batch, ncol=2, scales = "free")

Remove wells B12 and C12 from experiment batch 1A1

blanks0 %>%
  mutate(expwell=paste(Batch, well, sep=".")) %>%
  filter(expwell!="1A1.B12", expwell!="1A1.C12") %>%
  ggplot(aes(x=Time, y=OD, color=well)) +
  geom_point() +
  facet_wrap(~Batch, ncol=2, scales = "free_x")

Assess blank variability

blankSummaryTime <- blanks0 %>%
  mutate(expwell=paste(Batch, well, sep=".")) %>%
  filter(expwell!="1A1.B12", 
         expwell!="1A1.C12") %>%
  with_groups(c(Time, Batch), summarize, timeRange=max(OD)-min(OD))

blankSummaryWell <- blanks0 %>%
  mutate(expwell=paste(Batch, well, sep=".")) %>%
  filter(expwell!="1A1.B12", 
         expwell!="1A1.C12") %>%
  with_groups(expwell, summarize, wellRange=max(OD)-min(OD))

timeBlankPlot <- blankSummaryTime %>%
  ggplot(aes(x="Range per Timepoint", y=timeRange)) +
  geom_violin(alpha=0) +
  geom_quasirandom(alpha=0.5) +
  stat_summary(fun = "mean", geom="crossbar", width=0.8, color="red", linetype=2) +
  scale_y_continuous(breaks = seq(0, 0.0125, 0.0025), limits = c(0, 0.0125)) +
  labs(x=NULL, y="OD Range")
  
wellBlankPlot <- blankSummaryWell %>%
  ggplot(aes(x="Range per Well", y=wellRange)) +
  geom_violin(alpha=0) +
  geom_quasirandom(alpha=0.5) +
  scale_y_continuous(breaks = seq(0, 0.0125, 0.0025), limits = c(0, 0.0125)) +
  stat_summary(fun = "mean", geom="crossbar", width=0.8, color="red", linetype=2) +
  theme(axis.text.y = element_blank()) +
  labs(x=NULL, y=NULL)
  
plot_grid(timeBlankPlot, wellBlankPlot)

summary(blankSummaryTime$timeRange)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.000000 0.002000 0.003000 0.003039 0.004000 0.009000

Limit of detection appears to be about 0.003.

blanks <- curve_raw_data %>%
  mutate(expwell=paste(Batch, well, sep=".")) %>%
  filter(Number=="Blank",
                  expwell!="1A1.B12", 
                  expwell!="1A1.C12") %>%
  with_groups(c(Batch, Time), summarize, blankOD=mean(OD))

curve <- curve_raw_data %>%
  filter(Number!="Blank") %>%
  left_join(blanks, by=c("Time", "Batch")) %>%
  mutate(OD=OD-blankOD) %>%
  select(-blankOD) %>%
  mutate(Sample=paste(Strain, Batch, BioRep, sep="."))

Clean and filter data

curve %>%
  filter(Sample %!in% midLogFail) %>% #filtering based on mid log failure
  ggplot(aes(x=Time, y=OD, color=Species, shape=BioRep, linetype=BioRep)) +
  geom_point(alpha=0.5) +
  scale_x_continuous(breaks = seq(0,48,12)) +
  scale_y_log10() +
  facet_wrap(~Strain) +
  theme(axis.text = element_text(size=9),
        axis.title = element_text(size=11),
        strip.text = element_text(size=10),
        legend.title = element_text(size=11)) +
  labs(x="Time (Hours)", y=bquote(OD[600]), color="Species", shape = "Biological Replicate", linetype="Biological Replicate")

In addition to samples removed for quality control at mid log step:

1A1: Remove well E9 (one technical replicate for 14018 biological replicate B)

1B1: Remove JCP8108 biological replicate A

curveFilt <- curve %>%
  mutate(ExpWell=paste(Batch, well, sep=".")) %>%
  filter(ExpWell!="1A1.E9",
         Sample %!in% c(midLogFail, "JCP8108.1B1.A"))

curveMeans <-  curveFilt %>%
  filter(OD>0) %>% # remove any technical replicates below detection
  with_groups(c(Sample, Number, BioRep, Strain, Species, Time, Batch), summarize, mOD=mean(OD), sd=sd(OD)) %>%
  with_groups(c(Sample, Batch), nest, Time=Time, mOD=mOD, sd=sd)

Calculate growth rates

Use rolling regression to calculate growth curves

# create the rolling regression function
roll_regress <- function(x){
  temp <- data.frame(x)
  mod <- lm(temp)
  temp <- data.frame(slope = coef(mod)[[2]],
                     slope_lwr = confint(mod)[2, ][[1]],
                     slope_upr = confint(mod)[2, ][[2]],
                     intercept = coef(mod)[[1]],
                     rsq = summary(mod)$r.squared, stringsAsFactors = FALSE)
  return(temp)
}
rollRegressData <- curveFilt %>%
  mutate(SampWell=paste(Sample, well, sep=".")) %>%
  arrange(SampWell, Time) %>%
  mutate(Strain=factor(Strain, levels=strains),
         Species=factor(Species, levels=species))

Choose window parameters

First look at curves with smallest log phases

curveFilt %>%
  filter(Strain %in% c("14018", "315-A", "C0011E4", "C0093B3")) %>%
  ggplot(aes(x=Time, y=OD)) +
  geom_point(alpha=0.5) +
  facet_wrap(~Strain, scales = "free_x") +
  scale_y_log10()

Choose 2 hours for rolling window size. Choose 0.006 as minumum filtering threshold.

# define window - here every 2 hours with measurements every 0.5 hours
num_points <- ceiling(2*60/(60*0.5)) 
num_points
## [1] 4

Calculate growth rates

# run rolling regression on lnOD ~ time
models <- rollRegressData %>%
  filter(OD>0.006) %>%
  mutate(ln_od=log(OD)) %>%
  group_by(SampWell) %>%
  do(cbind(model = select(., ln_od, Time) %>% 
           zoo::rollapplyr(width = num_points, roll_regress, by.column = FALSE, fill = NA, align = 'center'),
           Time = select(., Time),
           ln_od = select(., ln_od))) %>%
  rename_all(., gsub, pattern = 'model.', replacement = '')

# create predictions
preds <- models %>%
  filter(., !is.na(slope)) %>%
  group_by(SampWell, Time) %>%
  do(data.frame(time2 = c(.$Time - 2, .$Time + 2))) %>%
  left_join(., models) %>%
  mutate(pred = (slope*time2) + intercept)

# get growth rates
growth_rate <- models %>%
  group_by(SampWell) %>%
  filter(slope == max(slope, na.rm = TRUE))
rollRegressPlotSubset <- rollRegressData %>%
  group_by(Strain) %>%
  sample_n(1, replace = FALSE) %>%
  mutate(Strain=factor(Strain, levels=strains)) %>%
  arrange(Strain) %>%
  .$SampWell %>%
  unique

# plot rolling regression
rollRegressData %>%
  filter(OD>0.006,
         SampWell %in% rollRegressPlotSubset) %>%
  mutate(ln_od=log(OD),
         Strain=factor(Strain, levels=strains),
         SampWell=factor(SampWell, levels=rollRegressPlotSubset)) %>%
ggplot(aes(x=Time, y=ln_od)) +
  geom_point(size=1) +
  geom_line(aes(time2, pred, group = Time), col = 'red', (subset(preds, SampWell %in% rollRegressPlotSubset) %>% mutate(SampWell=factor(SampWell, levels=rollRegressPlotSubset)))) +
  geom_segment(aes(x = Time, y = ln_od-3, xend = Time, yend = ln_od), (subset(growth_rate, SampWell %in% rollRegressPlotSubset) %>% mutate(SampWell=factor(SampWell, levels=rollRegressPlotSubset)))) +
  geom_segment(aes(x = ln_od-3, y = ln_od, xend = Time, yend = ln_od), (subset(growth_rate, SampWell %in% rollRegressPlotSubset) %>% mutate(SampWell=factor(SampWell, levels=rollRegressPlotSubset)))) +
  geom_text(aes(label= Strain, x=16, y=1.5), size=3, hjust="center", check_overlap = TRUE) +
  facet_wrap(~SampWell) +
  ylim(-8.5, 2.5) +
  theme(strip.text = element_blank(),
        strip.background = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        panel.grid = element_blank(),
        panel.border = element_rect(color = "black", size=1)) +
  labs(x = NULL, y = NULL) -> growthCurveFigures
growthCurveFigures

Plot Growth Rates

ratePlot <- growth_rate %>%
  separate(SampWell, c("Strain", "Batch", "BioRep", "Well"), sep="\\.", remove=FALSE) %>%
  mutate(Sample=paste(Strain, Batch, BioRep, sep=".")) %>%
  left_join(unique(curve[,c("Strain", "Species")])) %>%
  with_groups(c(Sample, Species, Strain, Batch, BioRep), summarize, mSlope=mean(slope), sdSlope=sd(slope)) %>%
  mutate(Strain=factor(Strain, levels=strains),
         Species=factor(Species, levels=species)) %>%
  ungroup %>%
  ggplot(aes(x=Strain, y=mSlope)) +
  geom_point(alpha=0) +
  geom_boxplot(position = position_dodge2(width=1), alpha=0, color="black", show.legend = FALSE) +
  geom_pointrange(aes(x=Strain, y=mSlope, ymin=mSlope-sdSlope, ymax=mSlope+sdSlope, color=Species), size=0.5, alpha=0.6, position = position_dodge2(width=1), show.legend = FALSE) +
  ylim(0,1.5) +
  speciesColor +
  theme(axis.text.x = element_text(angle=45, hjust=1),
        axis.title.y = element_text(size=10),
        panel.grid = element_blank(),
        legend.position = "none") +
  labs(x=NULL, y="Growth Rate")


protocolBar0 <- midlog_df %>%
  select(Strain, Species, Protocol) %>%
  unique %>%
  ggplot() +
    geom_bar(aes(x=Strain, y=2, fill=Protocol), stat="identity", width = 1.1) +
    scale_fill_manual(values=c("cornflowerblue", "burlywood1", "forestgreen")) +
    theme_void() +
    theme(legend.position = "none")

speciesBar0 <- midlog_df %>%
  select(Strain, Species, Protocol) %>%
  unique %>%
  ggplot() +
    geom_bar(aes(x=Strain, y=1, fill=Species), stat="identity", width=1.1) +
    speciesFill +
    theme_void() +
    theme(legend.position = "none")

#Legends
protocol_legend <- cowplot::get_legend(
    protocolBar0 + 
    guides(color = guide_legend(nrow = 1)) +
    theme(legend.position = "bottom"))

species_legend <- cowplot::get_legend(
    speciesBar0 +
    guides(color=guide_legend(ncol = 1)) +
    theme(legend.position = "bottom"))

#plots
plots <- plot_grid(protocolBar0, speciesBar0, ratePlot, ncol = 1, rel_heights = c(0.08, 0.08, 1), align = "v", axis = "bt")
plot_grid(protocol_legend, species_legend, plots, ncol=1, rel_heights = c(0.1, 0.25,  1))

Compare Growth Rates

Isolate ANOVA

  1. Assumption of equal variances
  2. Assumption of normality
qqplot <- growth_rate %>%
  separate(SampWell, c("Strain", "Batch", "BioRep", "Well"), sep="\\.", remove=FALSE) %>%
  mutate(Sample=paste(Strain, Batch, BioRep, sep=".")) %>%
  left_join(unique(curve[,c("Strain", "Species")])) %>%
  with_groups(c(Sample, Species, Strain, Batch, BioRep), summarize, mSlope=mean(slope), sdSlope=sd(slope)) %>%
  ggplot(aes(sample=mSlope)) +
  geom_qq_line() +
  geom_qq()
qqplot

One-way ANOVA

strain_growth_rate_DF <- growth_rate %>%
  separate(SampWell, c("Strain", "Batch", "BioRep", "Well"), sep="\\.", remove=FALSE) %>%
  mutate(Sample=paste(Strain, Batch, BioRep, sep=".")) %>%
  left_join(unique(curve[,c("Strain", "Species")])) %>%
  with_groups(c(Sample, Species, Strain, Batch, BioRep), summarize, mSlope=mean(slope), sdSlope=sd(slope)) %>%
  mutate(Strain=factor(Strain, levels=strains))

strain_growth_rate_DF %>%
  anova_test(mSlope~Strain) %>%
  formattable()
Effect DFn DFd F p p<.05 ges
Strain 21 61 14.494 1.79e-16
0.833

Tukey’s post-hoc test

growthRatePostHoc <- strain_growth_rate_DF %>%
  tukey_hsd(mSlope~Strain, p.adjust.method = "BH")

# growthRatePostHoc %>%
#  arrange(p.adj) %>%
# write.xlsx(file=suppTableOut, sheetName="TableS3_growth_rate_posthoc", col.names=TRUE, append=FALSE)
min(growthRatePostHoc$p.adj)

pvals <- growthRatePostHoc %>%
  select(group1, group2, p.adj) %>%
  dplyr::rename(var1=group2, var2=group1, cor=p.adj) %>%
  cor_spread() %>%
  column_to_rownames(var="rowname") %>%
  as.matrix()

growthRatePostHoc %>%
  select(group1, group2, estimate) %>%
  dplyr::rename(var1=group2, var2=group1, cor=estimate) %>%
  cor_spread() %>%
  column_to_rownames(var="rowname") %>%
  as.matrix() %>%
  corrplot(is.corr=FALSE, type = "lower", tl.col="black", p.mat = pvals, insig = "label_sig", sig.level = c(0.001, 0.01, 0.05), pch.cex = 0.9)

#bg="light gray", cl.pos="n",
grid.echo()

growthRatePostHocPlot <- grid.grab()

Moran’s I

Load tree

We choose the weights as \(_{wij} = 1/d_{ij}\) , where the d’s is the distances measured on the tree:

We can now perform the analysis with Moran’s I:

Carrying Capacity

carryingCapacity <- curveFilt %>%
  with_groups(c("well", "Number", "BioRep", "Species", "Strain",  "Batch",   "Sample",  "ExpWell"), summarize, OD=max(OD)) %>%
  with_groups(c("Number", "BioRep", "Species", "Strain",  "Batch",   "Sample"), summarize, maxOD=mean(OD), maxODsd=sd(OD)) %>%
  mutate(Strain=factor(Strain, levels=strains),
         Species=factor(Species, levels=species))

Plot

carryingCapacityPlot <- carryingCapacity %>%
  mutate(Strain=factor(Strain, levels=strains)) %>%
  ggplot(aes(x=Strain, y=maxOD)) +
  geom_boxplot(alpha=0) +
  geom_pointrange(aes(ymin=maxOD-maxODsd, ymax=maxOD+maxODsd, color=Species), size=0.5, alpha=0.6, position = position_dodge2(width=1), show.legend = FALSE) +
  speciesColor +
  theme(axis.text.x = element_text(angle=45, hjust=1),
        axis.title.y = element_text(size=10),
        panel.grid = element_blank(),
        legend.position = "none") +
  labs(x=NULL, y=bquote("Carrying Capacity (OD"[600]~")"))

plots2 <- plot_grid(protocolBar0, speciesBar0, carryingCapacityPlot, ncol = 1, rel_heights = c(0.08, 0.08, 1), align = "v", axis = "bt")
plot_grid(protocol_legend, species_legend, plots2, ncol=1, rel_heights = c(0.1, 0.25,  1))

Compare Carrying Capacity

Isolate ANOVA

  1. Assumption of equal variances
  2. Assumption of normality
qqplot <- carryingCapacity %>%
  ggplot(aes(sample=maxOD)) +
  geom_qq_line() +
  geom_qq()
qqplot

One-way ANOVA

strain_growth_rate_DF <- carryingCapacity %>%
  mutate(Strain=factor(Strain, levels=strains))

strain_growth_rate_DF %>%
  anova_test(maxOD~Strain) %>%
  formattable()
Effect DFn DFd F p p<.05 ges
Strain 21 61 6.81 1.96e-09
0.701

Tukey’s post-hoc test

carCapacPostHoc <- carryingCapacity %>%
  tukey_hsd(maxOD~Strain, p.adjust.method = "BH")

# carCapacPostHoc %>%
#  arrange(p.adj) %>%
#  write.xlsx(file=suppTableOut, sheetName="TableS4_car_capac_posthoc", col.names=TRUE, append=TRUE)
min(carCapacPostHoc$p.adj)

pvals <- carCapacPostHoc %>%
  select(group1, group2, p.adj) %>%
  dplyr::rename(var1=group2, var2=group1, cor=p.adj) %>%
  cor_spread() %>%
  column_to_rownames(var="rowname") %>%
  as.matrix()

carCapacPostHoc %>%
  select(group1, group2, estimate) %>%
  dplyr::rename(var1=group2, var2=group1, cor=estimate) %>%
  cor_spread() %>%
  column_to_rownames(var="rowname") %>%
  as.matrix() %>%
  corrplot(is.corr=FALSE, type = "lower", tl.col="black", p.mat = pvals, insig = "label_sig", sig.level = c(0.001, 0.01, 0.05), pch.cex = 0.9)

grid.echo()

carCapacPostHocPlot <- grid.grab()

Moran’s I

We can now perform the analysis with Moran’s I:

Figure 3: Growth stats

figure3B <- plot_grid(protocolBar0, speciesBar0, ratePlot + theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()), carryingCapacityPlot, ncol = 1, rel_heights = c(0.08, 0.08, .75, 1), align = "v", axis = "bt")
plot_grid(ggplot()+geom_blank()+theme_void(), growthCurveFigures, figure3B, rel_widths = c(0.08, 1.2, 1), nrow = 1, labels = c("A", "", "B"), label_size = 15)

#ggsave(file.path(figureOut, paste(Sys.Date(), "Figure3_growth_rates.png", sep = "_")))

Figure S2: Post Hoc Test Results

plot_grid(ggplot()+geom_blank(), growthRatePostHocPlot, ggplot()+geom_blank(), carCapacPostHocPlot, nrow = 1, rel_widths = c(0.1, 1, 0.1, 1), labels = c(" ", "A", " ", "B"), label_size = 20)

#ggsave(file.path(figureOut, paste(Sys.Date(), "FigureS2_growth_posthocs.png", sep = "_")))

CFUs

cfuTable1A1 <- tribble(
  ~Number, ~BioRep, ~Dilution, ~Colonies,
  "1", "A", 1e-4, 34, 
  "1", "B", 1e-4, 132,
  "2", "A", 1e-4, 70,
  "2", "B", 1e-3, 166,
  "3", "A", 1e-3, 83,
  "3", "B", 1e-2, 242,
  "4", "A", 1e-4, 32,
  "4", "B", 1e-3, 31,
  "5", "A", 1e-4, 48, 
  "5", "B", 1e-4, 23,
  "6", "A", 1e-4, 76, 
  "6", "B", NA, NA,
  "7", "A", 1e-3, 110,
  "7", "B", 1e-4, 71,
  "8", "A", 1e-3, 76,
  "8", "B", 1e-5, 25)

cfuTable1A2 <- tribble(
  ~Number, ~BioRep, ~Dilution, ~Colonies,
  "1", "A", 1e-3, 53, 
  "1", "B", 1e-4, 51,
  "2", "A", 1e-3, 29,
  "2", "B", 1e-3, 103,
  "3", "A", 1e1, 0,
  "3", "B", 1e1, 0,
  "4", "A", 1e-3, 115,
  "4", "B", 1e-3, 31,
  "5", "A", 1e-3, 127, 
  "5", "B", 1e-4, 26,
  "6", "A", 1e-3, 106, 
  "6", "B", 1e-2, 132,
  "7", "A", 1e-3, 47,
  "7", "B", 1e-2, 76,
  "8", "A", 1e-3, 60,
  "8", "B", 1e-3, 42)

cfuTable1B1 <- tribble(
  ~Number, ~BioRep, ~Dilution, ~Colonies,
  "1", "A", 1e-4, 5, 
  "1", "B", 1e-5, 2,
  "2", "A", 1e-1, 69,
  "2", "B", 1e-1, 74,
  "3", "A", 1e-4, 36,
  "3", "B", 1e-4, 35,
  "4", "A", 1e-3, 25,
  "4", "B", 1e-4, 28,
  "5", "A", 1e-3, 10, 
  "5", "B", 1e-4, 21,
  "6", "A", 1e-3, 32, 
  "6", "B", 1e-3, 26,
  "7", "A", 1e-3, 44,
  "7", "B", 1e-3, 30,
  "8", "A", 1e-4, 27,
  "8", "B", 1e-4, 24)

cfuTable1B2 <- tribble(
  ~Number, ~BioRep, ~Dilution, ~Colonies,
  "1", "A", 1e-3, 20, 
  "1", "B", 1e-3, 8,
  "2", "A", 1e-3, 21,
  "2", "B", 1e-3, 3,
  "3", "A", 1e-4, 58,
  "3", "B", 1e-4, 26,
  "4", "A", 1e-4, 25,
  "4", "B", 1e-2, 27,
  "5", "A", 1e-3, 59, 
  "5", "B", 1e-3, 14,
  "6", "A", 1e-3, 53, 
  "6", "B", 1e-3, 36,
  "7", "A", 1e-3, 37,
  "7", "B", 1e-3, 21,
  "8", "A", 1e-1, 84,
  "8", "B", 1e-3, 80,
  "9", "A", 1e-2, 6,
  "9", "B", 1e-1, 69)

cfuTable2A1 <- tribble(
  ~Number, ~BioRep, ~Dilution, ~Colonies,
  "1", "A", 1e-4, 26, 
  "1", "B", 1e-4, 87,
  "2", "A", 1e-4, 45,
  "2", "B", 1e-5, 39,
  "3", "A", NA, NA,
  "3", "B", NA, NA)

cfuTable2A2 <- tribble(
  ~Number, ~BioRep, ~Dilution, ~Colonies,
  "1", "A", 1e-5, 22, 
  "1", "B", 1e-5, 17,
  "2", "A", 1e-4, 6,
  "2", "B", 1e1, 0,
  "3", "A", NA, NA,
  "3", "B", NA, NA)

cfuTable3A1 <- tribble(
  ~Number, ~BioRep, ~Dilution, ~Colonies,
  "1", "A", 1e-2, 66, 
  "1", "B", 1e-3, 12,
  "2", "A", 1e-2, 73,
  "2", "B", 1e-2, 95,
  "3", "A", 1e-3, 33,
  "3", "B", 1e-3, 22,
  "4", "A", NA, NA,
  "4", "B", NA, NA)

cfuTable3A2 <- tribble(
  ~Number, ~BioRep, ~Dilution, ~Colonies,
  "1", "A", 1e-2, 66, 
  "1", "B", 1e-3, 12,
  "2", "A", 1e-2, 73,
  "2", "B", 1e-2, 95,
  "3", "A", 1e-3, 33,
  "3", "B", 1e-3, 22,
  "4", "A", NA, NA,
  "4", "B", NA, NA)


cfuTables <- list(cfuTable1A1, cfuTable1A2, cfuTable1B1, cfuTable1B2, cfuTable2A1, cfuTable2A2, cfuTable3A1, cfuTable3A2) %>%
  map(~mutate(.x, CFU_mL=(Colonies)/(Dilution*0.007))) %>%
  map2(., batches, ~mutate(.x, Batch=.y)) %>%
  purrr::reduce(full_join, by = c("Number", "BioRep", "Dilution", "Colonies", "CFU_mL", "Batch"))
strainKeys <- midlog_raw_data %>% # get strain information
  select(Number, Batch, Strain, Species) %>%
  unique()

cfuTables %>%
  left_join(strainKeys) %>%
  mutate(Protocol=str_sub(Batch, 1, 1),
         BatchNum=str_sub(Batch, 3,3),
         logCFU=log10(CFU_mL),
         Strain=factor(Strain, levels=strains),
         Species=factor(Species, levels=species)) %>%
  ggplot(aes(x=Strain, y=logCFU, color=Species, fill=Species, shape=BatchNum, linetype=BatchNum)) +
  geom_point() +
  stat_summary_bin(fun = "mean", geom = "crossbar", show.legend = FALSE) +
  scale_y_continuous(breaks = seq(1,10,1), limits = c(1,10), labels=c("10", expression(10^2), expression(10^3), expression(10^4), expression(10^5), expression(10^6), expression(10^7), expression(10^8), expression(10^9), expression(10^10))) +
  theme(axis.text.x = element_text(angle = 45, hjust=1),
        panel.grid.minor.y = element_blank(),
        legend.position = "bottom") +
  guides(shape=FALSE) +
  labs(x=NULL, y="CFU/mL")

Session Info

sessionInfo()
## R version 4.4.0 (2024-04-24)
## Platform: x86_64-apple-darwin20
## Running under: macOS Ventura 13.6.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] gridGraphics_0.5-1 ggpubr_0.6.0       cowplot_1.1.3      ggbeeswarm_0.7.2  
##  [5] corrplot_0.92      formattable_0.2.1  ape_5.8            rstatix_0.7.2     
##  [9] broom_1.0.6        xlsx_0.6.5         lubridate_1.9.3    forcats_1.0.0     
## [13] stringr_1.5.1      dplyr_1.1.4        purrr_1.0.2        readr_2.1.5       
## [17] tidyr_1.3.1        tibble_3.2.1       ggplot2_3.5.1      tidyverse_2.0.0   
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.5      beeswarm_0.4.0    xfun_0.46         bslib_0.7.0      
##  [5] htmlwidgets_1.6.4 rJava_1.0-11      lattice_0.22-6    tzdb_0.4.0       
##  [9] vctrs_0.6.5       tools_4.4.0       generics_0.1.3    parallel_4.4.0   
## [13] fansi_1.0.6       highr_0.11        pkgconfig_2.0.3   lifecycle_1.0.4  
## [17] farver_2.1.2      compiler_4.4.0    munsell_0.5.1     carData_3.0-5    
## [21] vipor_0.4.7       htmltools_0.5.8.1 sass_0.4.9        yaml_2.3.9       
## [25] crayon_1.5.3      pillar_1.9.0      car_3.1-2         jquerylib_0.1.4  
## [29] cachem_1.1.0      abind_1.4-5       nlme_3.1-165      tidyselect_1.2.1 
## [33] digest_0.6.36     stringi_1.8.4     labeling_0.4.3    fastmap_1.2.0    
## [37] colorspace_2.1-0  cli_3.6.3         magrittr_2.0.3    xlsxjars_0.6.1   
## [41] utf8_1.2.4        withr_3.0.0       scales_1.3.0      backports_1.5.0  
## [45] bit64_4.0.5       timechange_0.3.0  rmarkdown_2.27    bit_4.0.5        
## [49] ggsignif_0.6.4    zoo_1.8-12        hms_1.1.3         evaluate_0.24.0  
## [53] knitr_1.48        rlang_1.1.4       Rcpp_1.0.13       glue_1.7.0       
## [57] vroom_1.6.5       rstudioapi_0.16.0 jsonlite_1.8.8    R6_2.5.1